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Abstract 

The non-linear Poisson-Boltzmann (PB) equation for a circular, uniformly 
charged platelet, confined together with co- and counter-ions to a cylindrical 
cell, is solved semi-analytically by transforming it into an integral equation 
and solving the latter iteratively. This method proves efficient and robust, and 
can be readily generalized to other problems based on cell models, treated 
within non-linear Poisson-like theory. The solution to the PB equation is 
computed over a wide range of physical conditions, and the resulting osmotic 
equation of state is shown to be in fair agreement with recent experimental 
data for Laponite clay suspensions, in the concentrated gel phase. 
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Upon the addition of water, initially dry clays swell into stacks of parallel charged 
platelets separated by increasingly thick layers of water containing microscopic co- and 
counter-ions. This simple picture becomes less realistic as the concentration by weight of 
clay decreases and the initial smectic orientational order is gradually lost in favor of a dis- 
ordered gel structure. The development of this gel phase is responsible for characteristic 
rheological properties of the clay dispersion, which are very sensitive to salt concentra- 
tion and of considerable technological importance, in particular for the oil drilling indus- 
try. Recent experimental work (!]-[§ focussed on correlating the rheological behavior with 
the mesoscopic structure of the well-characterized synthetic clay Laponite, made of quasi- 
monodisperse disc-like platelets of thickness ~ lnm, diameter a 30 nm, and carrying a 
negative surface charge — Ze « — 10 3 e (where e is the proton charge). The complex interplay 
of excluded volume and electrostatic interactions between such highly anisotropic particles 
makes a theoretical description of clay gels extremely difficult [|J . However, the early stages 
of swelling are described reasonably well by a cell model where each platelet, assumed to be 
an infinite uniformly charged plane, is confined with its monovalent co- and counter-ions to 
a slab, and the distribution of the latter is treated within Poisson-Boltzmann (PB) theory 
0. As swelling proceeds the spacing between platelets becomes a sizeable fraction of their 
lateral dimension. Consequently, edge (finite size) effects become important and the prob- 
lem ceases to be one-dimensional. In the extreme case of an isolated platelet immersed in 
an electrolyte solution (infinite dilution limit), the PB problem can be solved numerically 
H|7]]. So far, these results provide the only available data for disc-like clays and no solution 
could be obtained for finite concentrations of platelets. In this Letter we report results for 
a cylindrical cell model, which accounts for both the finite size of the platelets and their 
finite concentration, and allows a semi-analytical treatment within non-linear PB theory. 
The linearized version of PB theory (LPB) for a circular platelet confined with its co- and 
counter-ions to the same cell has been recently considered in detail ||; however, LPB the- 
ory is not suited for large surface potentials, typically for Z > 10 2 , which is an order of 
magnitude down from physical Laponite charges, where linearization is no longer justified. 
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Motivated by the results of ||, we introduce the Wigner-Seitz cell around a clay platelet, 
and we consider that the average shape of this confining cell is cylindrical in the gel phase. 
If n is the number concentration of platelets, the volume of the cylindrical cell is v — 1/n — 
2nR 2 h where R is the radius and 2h the height, which is the spacing between platelets in the 
stack (typically of the order of 10 nm; however h is not a priori prescribed in our model). 
The platelets are assumed to be infinitely thin, of radius ro, and carrying a uniform surface 
charge density a = —Zejur^. Let p + (r) and p~(r) denote the counter-ion and co-ion density 
profiles inside the cell. The electrostatic potential <p(r) satisfies Poisson's equation: 



where <?p(r) = o~5(z)9(r — r) is the charge density of the platelet in cylindrical coordinates, 
and e is the macroscopic dielectric constant of water (primitive model). Within mean-field, 
or PB theory, the density profiles are approximated by the Boltzmann factors: 



where (3 = 1/ksT is the inverse temperature. The pre-factors p have individually no 
physical significance, while the product p$ Po depends on the conditions under consideration. 
If the salt concentration ns is assumed to be fixed in the dispersion, then the concentration of 
co- and counter-ions in the cell are imposed: n„ = N_/v = ns and n + = N + /v = n$ + Z/v, 
and the p$ are determined by the normalization constraints (canonical description): 



On the other hand, if the dispersion is in equilibrium with an ionic solution of concentration 
n' s , which acts as a reservoir, then the sum of the electrochemical potentials of co- and 
counter-ions in the cell must equal the same quantity for the salt in the reservoir. Since, 
for consistency, the ions in the cell are assumed to behave as an ideal solution, chemical 
equilibrium implies that H (semi-grand-canonical description): 




(1) 




(2) 




(3) 




(4) 
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By substituting (|2]) into (|T|), and subtracting K 2 ip(r) (where K is an arbitrary inverse 
length), we arrive at the closed non-linear partial differential equation for the potential 
<p{r): 

(V 2 - k 2 ) <p(r) = (g P (r) + ep+e-^) - e^e^W + ^<p(r)) (5) 

Bearing upon the physical significance of the cell model, it is assumed that the normal 
component of the electric field vanishes everywhere on the surface £ bounding the cylinder 
cell. Therefore, the boundary condition for the potential <p(r) is: 

[n(r)-V] reE ^(r)=0 (6) 

where n(r) denotes the normal to the surface £ at point r. The solution of (|5|) satisfies: 

<p(r) = ~\ G K (r, r') ^ P (r') + ep+e"^ - ep^e^ + ^(r')) & (?) 

where G K (r, r') is the Green's function which solves the linear problem: 

(V 2 -k 2 ) G K (r,r') = 5(r-r') (8) 

subject to the same boundary condition (Q). The reason for replacing the bare Laplace 
operator by the screened Laplace operator is that for k = the boundary condition (^) 
cannot be satisfied by the solution of eq. @, as becomes clear by integrating both sides of 
eq. (§) (with k = 0) over the cell volume v. The Neuman-like boundary condition allows 
the Green's function to be expanded, in cylindrical coordinates, in the form of a Bessel-Dini 
series, along the lines of the procedure for obtaining the solutions to the LPB problems in 
§. We find: 

G K (r,r') = £C(</>y)cosh 

n>l 

where r = (p,<ft,z), the + and — signs correspond to the situations z > z' and z < z' 
respectively, y n is the n th root of Ji(y) = 0, Jo and J\ are the zeroth and first order Bessel 
functions, A~ 2 = (y 2 + k 2 R 2 )/R 2 , and: 
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The kernel G K being explicitly known, the non-linear integral equation (|7|) may be solved for 
(p(r) by an iterative Picard procedure starting from an initial guess of f(r) (in most cases 
we have chosen either </?(r) = or the solution of LPB theory M). The arbitrary inverse 
length k has generally been chosen in the range [k d /10, 10 k d ], k d being the inverse Debye 
length in the reservoir. Cylindrical symmetry implies <£>(r) = (p(p,z), which is calculated 
on a two-dimensional n p x n z grid spanning half the cylinder, with n p = 200 grid points to 
cover the interval [0, R] (this must be sufficiently large to allow for a proper representation 
of the Bessel functions), and n z = 40 points for the interval [0, h]. 50 terms were retained 
in the Bessel-Dini series ([J), and at each step of the iteration typically 10% to 5% of the 
new estimate of ip(r) were mixed with 90% to 95% of the previous estimate to ensure proper 
convergence; the latter was generally achieved in about 50 to 200 iterations to ensure a 
relative accuracy of at least 10~ 5 . Examples of the resulting density profiles are shown in 
Fig. 1 and compared with the LPB predictions; under realistic conditions the deviations of 
LPB from PB theory are indeed considerable. Provided the iterative procedure converges, 
the method yields a solution independent of k, i.e. no limit k —>■ is required at the end. 
We have successfully implemented two independent checks for the accuracy of the computed 
potential. First, eq. (0) must be satisfied for any values of n. Secondly, since the multipoles 
of the total charge distribution inside the cell can be transformed from volume integrals {i.e. 
their definition) into surface integrals (making use of Poisson's equation), these two formulae 
must yield the same result provided the potential <£>(r) is a solution of the PB problem. 

Once the potential ip(r) is known for different surface charge densities a, the Helmholtz 
free energy F can be calculated by a charging process, e.g. at constant Debye length 0. In 
the semi-grand-canonical calculations, the number of co- and counter-ions, AL and N + , are 
computed from the density profiles (for a given salt concentration n' s in the reservoir), and 
the grand potential Q is estimated from: 

= F-2N„k B T\n \n' s X 3 } (11) 



where A is the de Broglie thermal wavelength. 

The total quadrupole moment Q of the charge distribution in the cell (first non-vanishing 
multipole), the osmotic pressure II and the disjoining pressure ILj, can also be calculated 
from the density profiles as described in ref. ||. For a given clay concentration n, and hence 
cell volume v, all these quantities depend on the aspect ratio h/R (or equivalently h/r ) of 
the cell. The equilibrium stacking is determined by minimizing either the free energy or the 
grand potential with respect to h/r . This is illustrated in Fig. 2. As expected, h is found to 
be of the order of 10 nm. The minima from the PB calculations turn out to be consistently 
flatter than the LPB results, pointing to the possibility of large fluctuations of the stacking 
configurations. The optimum aspect ratios h/r predicted by PB and LPB are practically 
identical. This ratio varies with n, but is very insensitive to salt concentration. As in the 
LPB case, Q is found to vanish at the minimum of the relevant thermodynamic potential 
(F in the canonical case or Q in a semi-grand-canonical description, the latter being suited 
for a comparison with the experiments of 0), while II and 11^ are equal at that point. This 
property is exact and may serve as a check for the consistency of the numerical calculations, 
while the vanishing of Q cannot be proven a priori within the PB framework, but may be 
regarded as an indication of the validity of our cell model for the description of dispersions 
of clay platelets ||. 

Our PB results for the osmotic equation of state, II versus clay concentration, are com- 
pared in Fig. 3 to the experimental data of Mourchid et al. for three reservoir salt 
concentrations n' s . The agreement is reasonable at the higher concentrations (above C ~ 
4% (w/w) in units of % of weight of Laponite per weight of solvent (water)) corresponding 
to the gel phase, except for the lowest salinity investigated (n' s = 10 _4 M with 1M= lmol 
dm -3 ) for which the convergence criterion is more difficult to meet. At this low salinity 
the PB predictions overestimate the osmotic pressure considerably, compared to the exper- 
iments. This disagreement may be linked to an excessive depletion of salt (Donnan effect) 
predicted by PB theory (see Table 1). Too strong a depletion leads to insufficient screening, 
and hence to an overestimation of the osmotic pressure. 
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The cell model should not apply at low clay concentrations where the parallel stacking 
of platelets is completely lost, giving rise to the orient at ionally disordered sol phase. Since 
there is no long-ranged orientational order in the gel phase our cell model tends to over- 
estimate the effect of lamellar order. However, the analytical results obtained in || within 
LPB indicate that given a clay concentration, and at least at this level of approximation, 
the precise shape of the confining cell does not affect significantly mesoscopic or macroscopic 
quantities such as the osmotic pressure. Finally, note that LPB theory leads to negative 
osmotic pressures under most conditions represented in the figure. 

In order to assess the importance of edge effects we have also considered the case of 
infinite planar platelets having the same charge and mass density than the finite Laponite 
discs. The corresponding one-dimensional PB problem can be solved numerically along the 
lines in ref. ||. For a given concentration by weight in Laponite the spacing 2h 00 between 
the infinite platelets differs from the optimum (lowest free energy) spacing between the finite 
discs (which is practically independent of n' s ). A comparison is made in Fig. 3. It is clear 
that as n' s increases edge effects become more important and the simple model of a stack 
of charged layers cannot describe the structure of the Laponite suspension. The two curves 
for the osmotic pressure as obtained from this simple model at the higher concentrations n' s 
are completely off scale. In view of the poor performance of PB theory for finite platelets 
at n' s = 10~ 4 M, the reasonable agreement of the osmotic pressure calculated for infinite 
platelets with experimental data must be considered as fortuitous. 

In summary, the present PB calculations confirm the qualitative features of earlier LPB 
results but while the latter predictions for the osmotic pressure are unphysical, the PB 
results are in reasonable agreement with experimental data for Laponite in the high concen- 
tration regime. The transformation of the initial Poisson equation into an integral equation 
problem, combined with the analytical calculation of the corresponding Green's function, 
proves efficient and can be readily generalized to other Poisson problems with non-linear 
source terms. 
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TABLES 



C % (w/w) 


n' s = 10~ 4 M 


Values of ns/n' s at 
n' s = 1Q- 3 M n' s = 5x 10~ 3 M 


n' s = 10- 2 M 


3.293 


0.025 


0.238 


0.654 


0.786 


4.390 


0.018 


0.168 


0.560 


0.716 


5.488 


0.013 


0.125 


0.476 


0.653 



TABLE I. Values for the ratio between the salt concentration in the system, n#, and in the 
reservoir, n's, for several Laponite concentrations. These are obtained from semi-grand canonical 
non-linear PB calculations, at the aspect ratio that minimizes the grand potential. 
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FIG. 1. Dimensionless electrostatic potential &(z, p) = (3eip(z,p) at z = versus p/£, and at 
p = versus z/l, with l~ 2 = 4-7T f3e 2 (2ns) / t-of-r ■ The profiles were obtained in non-linear PB, with 
the calculations performed in the canonical ensemble, and in LPB (see [|| for details). The results 
in this figure are for an aqueous solution (e r = 78) of clay discs of radius tq = 125 A and surface 
charge Z = 700, at a temperature T = 300 K. The concentration of clay is n = 5 x 10~ 5 M and that 
of added monovalent salt ns = 10~ 3 M. The aspect ratio is h/ro = 0.971, the value that minimizes 
the Helmholtz free energy (see Fig 2). 
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FIG. 2. a) Dependence of the dimensionless Helmholtz free energy (3(F — Fq)/Z and the di- 
mensionless quadrupole moment Q/Q dlsc , with Q dtsc = Zer^/A and Q = Q zz = —2Q XX = —2Q yy , 
on the aspect ratio h/ro of the Wigner-Seitz cylindrical cell. The solid lines represent (3(F — Fq)/Z 
and the dashed lines are for Q/Q dtsc . b) Dependence of the disjoining II^ and osmotic II pressure 
on h/ro- The dashed lines correspond to ILi and the solid lines to II, with IIo = kBT(2ns)- In 
both plots the curves with O were obtained from PB calculations whereas those with no symbols 
represent LPB results. At the optimum h/ro, which minimizes F, Q vanishes and II = ILj holds. 
The optimum value is consistently found to be h/ro = 0.971, and Fo is the value of F at this aspect 
ratio. The canonical ensemble calculations correspond to the same conditions as in Fig 1. The 
grand potential calculated in the semi-grand-canonical ensemble at a reservoir salt concentration 
n' s = 2.54257 x 10~ 3 M exhibits a minimum at the same value of h/ro, where Q vanishes and 



n = = 2530 Pa from the calculations in both ensembles. 
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FIG. 3. Equations of state of Laponite suspensions for different reservoir salt concentration n' s . 
The filled symbols represent experimental results obtained by Mourchid et al. 0. The solid lines 
with open symbols represent our semi-grand canonical PB results obtained at the same values of 
n' s . Each of these points corresponds to the osmotic pressure at the minimum of the grand potential 
with respect to h/r Q , at fixed n' s . The values of n' s are 10 _4 M (□), 1CT 3 M (y), 5 x 1(T 3 M (O) , 
and 10 _2 M (A). The dash-dotted, the dashed and dotted lines are the equations of state obtained 
from the one-dimensional PB calculation for a stack of infinite planar charged layers in the presence 
of salt, with the same reservoir concentration n' s , respectively. The computed results are for an 
aqueous solution (e r = 78) of Laponite at T = 300 K. The radius of the clay discs is ro = 150 A 
and the surface charge Z = 1000. 
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